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Acoustic  plane  wave  Grazing  angle  Reflection  coefficient 

Density  profile  Local  wave  Sea  bed 

Forward  scattering  Inverse  problem  Sound  speed  profile 

Geoacoustic  model _ Inversion  algorithm _ 
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--^Investigations  of  acoustic  bottom  reflectivity  attempt  to  infer  the 
structure  of  the  sea  bed  (e.g.,  the  density  and  sound  speed  profiles)  from  a 
limited  knowledge  of  the  reflection  coefficient.  For  many  applications,  an 
adequate  model  to  study  the  acoustic  interaction  is  provided  by  the  scatter¬ 
ing  of  plane  waves  from  a  one-dimensional  inhomogeneous  medium. 

In  contrast  to  formally  exact  solutions  to  this  inverse  scattering 
problem,  Candel  et  al .  (Journal  of  Sound  and  Vibration,  vol .  68,  1980, _ 

00,^1473 


20.  (Cont'd) 


pp.  571-595)  propose  an  approximate  scheme  that  can  be  readily  implemented 
numerically.  Their  method  applies  the  forward  scattering  approximation  to  a 
local  wave  decomposition  of  the  acoustic  field.  As  a  result,  the  plane  wave 
reflection  coefficient  is  obtained  as  a  nonlinear  Fourier  transform  of  the 
logarithmic  derivative  of  the  local  admittance.  Inversion  of  the  integral 
transform  enables  the  recovery  of  admittance  versus  depth  by  means  of  a  numer 
ical  integration  using  a  single  impulse  response  from  the  sea  bed.  Separate 
recovery  of  both  the  density  and  sound  speed  profiles  requires  at  least  two 
impulse  responses  corresponding  to  two  distinct  grazing  angles. 

This  paper  appraises  an  implementation  of  Candel  et  al.'s  inversion 
algorithm  for  the  recovery  of  the  density  and  sound  speed  profiles  from  two 
realistic  geoacoustic  models  of  the  sea  bed  for  which  bandlimited  impulse 
responses  were  synthetically  generated. 
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THE  DETERMINATION  OF  MATERIAL  PROPERTIES  OF  THE  SEA  BED 
FROM  THE  ACOUSTIC  PLANE  WAVE  REFLECTION  RESPONSE 

INTRODUCTION 


The  detection  capability  of  passive  sonar  systems  that  use  bottom 
bounce  propagation  paths  to  estimate  range  and  bearing  to  a  target  depends 
on  the  reflection  coefficient  of  the  sea  bed.  At  low  frequencies,  the 
reflection  process  can  be  complicated  by  long  wavelength  sound  waves  that 
penetrate  the  surflclal  sediments.  As  a  result,  proper  Interpretation  of 
the  acoustic  signals  following  subbottom  Interactions  requires  a  knowledge 
of  the  material  properties  (e.g.,  the  density  and  sound  speed  profiles) 
below  the  sea  bed.  A  basic  problem  of  underwater  acoustics  Is  how  to 
accurately  determine  the  structure  of  the  sea  bed  from  a  limited  knowledge 
of  Its  reflection  coefficient.  For  many  applications,  an  adequate  model  to 
study  this  Inverse  problem  Is  provided  by  the  scattering  of  acoustic  plane 
waves  from  a  one-dimensional  Inhomogeneous  medium. 

Formal  methods  for  solving  the  one-dimensional  Inverse  scattering 
problem  are  reviewed  by  Burrldge1  and  Newton.2  In  contrast  to  these 
exact  solutions  that  provide  few  numerical  results,  Candel,  DeFI 11 1  pi ,  and 
Launay2’4  propose  an  approximate  one  that  Is  readily  Implemented 
numerically.  In  their  approach,  the  acoustic  field  Is  decomposed  Into  the 
local  upgolng  and  downgoing  plane  waves  described  In  Claerbout.9 
Application  of  the  forward  scattering  approximation  leads  to  an  analytic 
representation  for  the  reflection  coefficient  as  a  nonlinear  Fourier 
transform  of  the  logarithmic  derivative  of  the  local  admittance.  Inversion 
of  the  Integral  transform  results  In  a  numerical  algorithm  by  which  the 
admittance  profile  of  the  Inhomogeneous  medium  can  be  recovered  from 
Integration  of  a  single  Impulse  response.  At  least  two  Impulse  responses 
corresponding  to  two  distinct  grazing  angles  of  Incidence  are  required  to 
recover  both  density  and  sound  speed  profiles. 

This  paper  reviews  the  basic  results  leading  to  the  Inversion  algorithm 
of  Candel  et  al.*»4  It  also  appraises  the  numerical  Implementation  of 
their  method  for  recovering  both  density  and  sound  speed  profiles  from  two 
realistic  geoacoustic  models  of  the  sea  bed  for  which  Impulse  responses  are 
synthetically  generated. 


BASIC  THEORY 

The  mathematical  model  Is  shown  In  figure  1.  The  region  0  <  z  <  H  Is 
occupied  by  an  Inhomogeneous  liquid  that  has  an  arbitrary  variation  of 
density,  p(z),  and  sound  speed,  c(z).  The  homogeneous  regions  are 
characterized  by  constant  density,  sound  speed  pairs:  pq,  c0  for  z 

<  0  and  px,  cx  for  z  >  H.  All  regions  are  assumed  to  be 

nonabsorbing. 

An  acoustic  plane  wave,  p^.  Initially  propagating  downward  at  an 
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Figure  1.  The  Mathematical  Model 


angle,  eQ,  to  the  z  ■  0  plane  encounters  the  Inhomogeneous  medium  at  z  * 
0.  Within  0  <  z  <  H,  the  Incident  wave  Is  partially  reflected,  and  a 

reflected  wave,  Rpr,  Is  returned  to  z  <  0.  In  z  >  H,  where  no  further 

reflections  occur,  only  a  transmitted  wave,  Tpt.  continues  to  propagate. 

For  stratified  media,  all  field  quantities  are  Independent  of  y  (say) 
and  exhibit  a  common  wavenumber  In  the  x-dlrectlon  (l.e.,  kx  * 
k0cose0).  This  number  Is  fixed  by  the  grazing  angle  of  the  Incident 

wave.  The  assumption  of  time-harmonic  waves  with  angular  frequency,  u  * 
2wf ,  allows  suppression  of  the  common  factor,  exp[1kx  x  -  1«t],  from 
all  field  quantities. 

The  analytical  route  taken  by  Candel  et  al.*»*  to  solve  the 
one-dimensional  Inverse  scattering  problem  Is  summarized  In  figure  2.  In 
step  (1),  the  plane  wave  reflection  coefficient,  R,  Is  obtained  formally  by 
numerical  Integration  of  the  basic  equations,  which  are  satisfied  by  the 

pressure,  p,  and  the  z-component  of  particle  velocity,  w.  Here.  Yz  Is  a 
"longitudinal  admittance,”  defined  as  Yz  *  (kz/k)/(pc)  «  kz/(<*>), 
where  kz  ■  [k*  -  kx*]1'*  Is  the  z-component  of  the  total 

wavenumber,  k  «  w/c,  at  a  given  depth  z.  The  splitting  Into  a  local 
upgolng  wave,  U,  and  a  downgoing  wave,  D,  shown  In  step  (2),  results  In  the 
local  wave  representation,  shown  In  step  (3).  Although  this  splitting  1$ 
not  unique,  as  described  by  Schelkunoff*  and  Sluljter,7  the  local 
reflection  coefficient,  defined  by  R  »  U/D,  gives  the  exact  result  for 
uniform  media.  For  |U|  «  |D|,  physical  considerations  suggest  that 

the  Interpretation  Is  suitable  for  Inhomogeneous  media.  Application  of  the 
forward  scattering  approximation  In  step  (4)  results  In  an  analytical 
representation  for  R  as  a  nonlinear  Fourier  transform  of  the  logarithmic 
derivative  of  the  local  admittance,  g.  For  |g|  «  1,  this  same  result 
was  obtained  by  Schelkunoff*  and  Brekhovsklkh.*  Transformation  to  the 
new  depth  coordinate,  z*.  In  step  (S)  results  In  a  standard  Fourier 
Integral.  Evaluation  of  the  time  response  of  R  In  step  (6)  determines  the 
basic  result  that  the  Impulse  response  at  time,  t.  Is  related  to  the 
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(1)  Global  field  representation  (2)  Local  wave  decomposition 


dp/dz  *  Iwpw,  p(H)»l 

dw/dz  -  1kzYzp,  w(H)-Yz(H) 

_ y 

U  ■  (p  -  w/Yz)/2 

p(0)  -  w(0)/Y  (0) 

R(0) . . 

P(0)  +  w(0)/Yz(0) 

0  -  (P  +  w/Yz)/2 

(4)  Forward  scattering  approximation 

! 

(3)  Local  wave  representation 

U<°>  -  0,  dO{0)/dz  -  [1k2  -  g]0^0) 

dD/dz  -  1kzD  -  g(D  -  U),  0(H)  «  1 

dU^/dz  »  -[1k2  +  g]U^  +  g 

dU/dz  -  -ikzU  ♦  g(D  -  U),  U(H)  «  0 

R(0)  «  U(1,(0)/D*0)(0) 

H  S 

■  - /g(s)exp[21/  k  (t)dt]ds 

0  0  z 

R(z)  •  U(z)/0(z) 

g  -  (dYz/dz)/(2Yz) 

! 

(5)  Standard  Fourier  transform 

(6)  Time  (Impulse)  response 

z 

2*  ■  (2/k  )  /k  (t)dt 

0  0  z 

r(t)  -  /R(f)exp(-2ir1ft)df 
-00 

Z*(H) 

R(f)  *  -  / g(z*)exp(2Tt1fz*/c0)dz* 

-  -c0(dYz/dz*)/(2Yz)| 

z*  -  cQt 

V 

(7)  Inversion  algorithm  for  adnlttance  vs  depth 

<iydz*  -  -(Z/c0}rY2.  Yz(0)  -  Slne0/p0c0 
dz/dz*  -  l/(2pc0Yz),  z(0)  -  0 


Figure  2.  Flowchart  of  the  Inversion  Algorithm 
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variation  In  material  properties  about  an  "active"  layer  at  depth,  z*  .« 
c0t.  Finally,  In  step  (7),  two  first  order  differential  equations  are 
obtained  for  the  numerical  determination  of  Yz  versus  z  from  a  single 
reflection  sequence,  r(t). 

Two  aspects  of  this  analysis  deserve  comment.  First,  the  manipulations 
Inherent  In  step  (6)  require  that  the  relative  variation  of  admittance  be 
Independent  of  frequency.  Although  this  Is  a  reasonable  assumption  for 
density  and  sound  speed,  usual  treatments  of  absorption  presume  a  linear 
dependence  on  frequency.  For  this  reason,  absorption  was  taken  to  be  zero 
In  the  mathematical  model.  Second,  although  the  forward  scattering 
approximation  leads  to  a  direct,  noniterative  Inversion  algorithm,  only 
primary  reflections  within  the  Inhomogeneous  medium  can  be  accommodated.  As 
a  result,  the  presence  of  subbottom  multiple  reflections  In  the  Impulse 
response  sequence  can  limit  the  depth  to  which  the  admittance  profile  can  be 
recovered  accurately. 

The  result  In  step  (7)  Indicates  that  more  Information  Is  required  to 
recover  both  density  and  sound  speed  profiles.  If  the  subscripts  i  and 
a  are  used  to  denote  quantities  corresponding  to  grazing  angles  6i 
and  e*.  where  e,  >  e, ,  then  It  can  be  shown  that  recovery  of 
both  p  and  n  =  cQ/c  can  be  obtained  from  the  numerical  Integration  of 
the  four  differential  equations. 


«i(Yz)a/<it*1  -  -<2/c0>  ri  (ya  . 

(1) 

dz/dz*a  -  [ZpCotY^r1  , 

(2) 

dz*a/dz*a  -  (Yz)a/(Yz)a  , 

(3) 

d(Yz),/dz*a  -  -<2/c„)  r,  (YZ)*/(YZ)X  , 

(*) 

together  with  the  relations 


A«*  -  [cosa  ea  -  cos*  ea]/[(Yz)a  -  (Yz)a  ]  ,  (5) 

n  -  C<V»  cos*  e>  _  cos*  ‘  (V*  ]  *  (6) 

Equations  (1)  through  (6)  are  the  basis  of  the  numerical  results  described 
In  the  remainder  of  this  paper.  Reflection  responses  are  required  for  two 
probing  directions,  and  four.  Instead  of  two,  differential  equations  with 
appropriate  Initial  conditions  must  be  Integrated. 


GEOACOUSTIC  MODELS 

The  density  and  sound  speed  profiles  of  the  two  geoacoustic  models  that 
provide  the  numerical  simulations  are  shown  In  figure  3.  Both  geoacoustic 
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models  were  determined  by  analysis  of  deconvolved  bottom  reflected  signals 
generated  by  explosive  signal,  underwater-sound  (SUS)  charges.  The 
measurements  for  the  Hatteras  Abyssal  Plain  (28*30'  N.,  70°30'  W.)  were  made 
by  the  Naval  Underwater  Systems  Center  (NUSC)  and  Interpreted  by  Hersteln, 
Dullea,  and  Santanlello.*  This  model  was  used  In  other  time-waveform 
simulations  by  OINapoll,  Potter,  and  Hersteln. The  measurements  for  the 
Tufts  Abyssal  Plain  (46*00'  N.,  143*30'  W.)  were  made  by  the  Oefence 
Research  Establishment  Pacific  (OREP)  and  Interpreted  by  Chapman.11  Both 
models  Incorporate  refracting  sediment  layers  separated  by  discontinuities 
In  material  properties. 


NUMERICAL  RESULTS 


The  numerical  examples  provided  by  Candel  et  al.4  were  restricted  to 
waves  propagating  at  normal  Incidence  (e0  *  90*)  In  regions  of  constant 
density  (p  ■  pe)  only.  For  their  results,  recovery  of  n(z)  required 
only  a  single  Impulse  response,  as  Indicated  In  step  (7)  of  figure  2. 
Thomson  and  Chow1*  extended  the  numerical  Implementation  to  Include  two 
Impulse  responses,  corresponding  to  normal  and  oblique  grazing  angles,  so 
that  simultaneous  recovery  of  both  p(z)  and  n(z)  Is  possible.  For  the 
geoacoustic  models  In  figure  3,  a  similar  approach  Is  required. 


For  each  geoacoustic  model,  the  computer  simulation  proceeded  as 
follows.  At  a  given  grazing  angle,  the  complex,  bandllmlted  frequency 
response,  R  ■  Re[R]  1Im[R],  was  computed  at  256  discrete  frequencies,  f^ 
-  kaf ,  k  ■  0,  1  ,  ....  255,  using  af  -  0.5  Hz.  The  sequence  R  was 

extended  to  N  ■  1024  points  by  appending  zero  values.  A  discrete  fast 
Fourier  transform  was  used  to  produce  1024  estimates  of  the  time  response, 

r,  at  t„  -  nat,  n  -  0,  1 .  1023,  where  at  *  l/(Naf)  =  1/512 

seconds.  Each  time  response  was  convolved  with  a  lowpass  digital  filter 
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Figure  4.  Frequency  Responses  of  the  Reflection  Coefficient, 
R  -  |R|exp[1arg(R)],  at  Grazing  Angles  of  60*  and  90* 
for  the  Hatteras  Abyssal  Plain  Model 


designed  using  the  window  method  as  described  by  Kaiser  and  Reed.1*  This 
filter  Inhibited  Gibbs  oscillations  caused  by  sampling  and  truncation 
effects.  For  the  results  In  this  paper,  a  31-point  Kaiser  window  with  a 
stopband  suppression  of  60  dB  and  a  cutoff  frequency  of  64  Hz  was  used. 

Each  filtered  time  response  was  multiplied  by  &f  to  approximate  the 
analytical  Fourier  transform  result. 

Figure  4  depicts  the  passband  frequency  responses  of  the  reflection 
coefficient  for  the  Hatteras  Abyssal  Plain  model.  The  amplitude  |R|, 
and  phase,  arg(R),  are  shown  for  grazing  angles  of  60*  and  90*.  |R| 

exhibits  two  distinct  modulations  caused  by  reflections  from  layers  of 

different  thicknesses.  Arg(R)  Increases  monotonlcally  but  is  nonlinear 
because  of  the  sound  speed  gradients  within  the  thick  layer.  The  greater 
amplitudes  that  occur  for  oblique  Incidence  than  for  normal  Incidence  are 
typical  for  Impedance  values  that  Increase  with  depth. 

The  filter*  Impuls  responses  for  this  model  are  shown  In  the  upper 

part  of  figure  b.  he  »ft  trace  corresponds  to  a  grazing  angle  of  60*,  and 
the  right  trace  Cvrresponds  to  a  grazing  angle  of  90*.  Each  response  Is 


Figure  5.  Filtered  Time  Responses  of  the  Reflection  Coefficients 
for  Grazing  Angles  of  60°  and  90°  (Top)  Used  to  Obtain  the 
Normalized  Oenslty  and  Sound  Speed  Reconstructions  for  the 
Hatteras  Abyssal  Plain  Model  (Bottom) 


normalized  by  Its  peak  value.  The  Initial  time  offset  for  each  trace 
corresponds  to  the  two-way  travel  time  between  z  *  0  and  a  receiver  located 
at  z  •  -100  m.  The  smaller  time  offset  observed  for  oblique  Incidence 
results  from  the  plane  wave  propagating  downward  at  the  longitudinal  phase 
speed,  «/kz.  The  three  "pulses"  observed  at  both  grazing  angles  are 
associated  with  discontinuities  at  the  layer  Interfaces,  z  ■  0,  20.4,  and 
300  m.  Between  the  second  and  third  reflections,  continuous  returns  are 
observed,  which  are  a  result  of  the  refracting  nonzero  sound  speed  gradients 
within  the  deep  layer. 

Use  of  these  filtered  time  responses  as  Input  to  equations  (1)  through 
(6)  determines  the  reconstructed  values  of  p/p0  and  n  «  c0/c,  shown 
In  the  lower  part  of  figure  5.  For  comparison,  the  density  and  sound  speed 
profiles  of  the  Hatteras  Abyssal  Plain  model  used  to  generate  the  synthetic 
data  are  displayed  on  the  same  graphs.  It  Is  evident  that  the  global 
behavior  of  the  reconstructed  profiles  agrees  with  the  model  Inputs.  For 
the  highest  passband  frequency  used  In  the  Inversion  process,  the  thin  upper 
layer  In  this  model  Is  only  0.8  wavelengths  thick.  At  this  wavelength,  the 
small  changes  In  n(z)  are  not  resolved,  although  some  of  the  difficulty  Is 
due  to  remanent  Gibbs  effects.  It  Is  Interesting  to  note  that  the  larger 
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Figure  6.  Frequency  Responses  of  the  Reflection  Coefficient, 
R  =  |R|exp[1arg(R)],  at  Grazing  Angles  of  70°  and  90“ 
for  the  Tufts  Abyssal  Plain  Model 


density  variations  within  this  layer  are  easily  resolved  by  the  bandlimlted 
time  responses. 

Figure  6  shows  the  frequency  responses  of  the  reflection  coefficients 
computed  for  the  Tufts  Abyssal  Plain  model.  The  large  Impedance  Increase  at 
the  base  of  the  Inhomogeneous  region  produces  values  of  |R|  greater  than 
those  found  for  the  Hatteras  Abyssal  Plain  model.  In  this  case,  the  forward 
scattering  criterion,  |U|  «  |D|  ,  Is  only  marginally  satisfied.  As 
before,  a  double  modulation  caused  by  interference  between  reflections  at 
the  top  and  bottom  of  each  Inhomogeneous  layer  occurs.  The  phase,  arg(R), 
which  changes  rapidly  with  increases  In  frequency,  Is  nonlinear  because  of 
the  refracting  sediments. 

The  filtered  time  responses,  corresponding  to  the  frequency  responses 
given  In  figure  6,  are  presented  In  the  upper  part  of  figure  7.  The  first 
three  "pulses"  correspond  to  primary  reflections  from  the  Interfaces,  z  «*  0, 
35,  and  200  m.  Multiple  subbottom  reflections  arriving  at  later  times  are 
clearly  visible.  The  lower  part  of  figure  7  shows  the  reconstructed  values 
of  density  and  sound  speed  compared  to  the  Input  values  for  the  Tufts 
Abyssal  Plain  model.  As  In  the  previous  example,  the  global  agreement 
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Figure  7.  Filtered  Time  Responses  of  the  Reflection  Coefficients 
for  Grazing  Angles  of  70*  and  90°  (Top)  Used  to  Obtain  the 
Normalized  Density  and  Sound  Speed  Reconstructions  for  the 
Tufts  Abyssal  Plain  Model  (Bottom) 


between  reconstructed  and  model  Input  profiles  Is  excellent.  Again, 
although  n(z)  within  the  acoustically  thin  upper  layer  Is  not  resolved  by 
limitations  of  the  64  Hz  band,  the  variations  In  density  of  this  region  are 
well  defined.  The  local  disagreement  near  z  =  200  m  apparently  results  from 
remanent  Gibbs  effects  caused  by  Imperfect  filtering.  For  this  model, 
reconstruction  of  the  density  and  sound  speed  profiles  to  greater  depths 
would  be  Inappropriate  because  of  multiple  subbottom  arrivals. 


SUMMARY 

This  document  has  reviewed  the  direct  Inversion  algorithm  proposed  by 
Candel  et  al.4  to  recover  the  density  and  sound  speed  profiles  of  the  sea 
bed  from  acoustic  reflection  data.  A  numerical  Implementation  was  applied 
to  the  reflection  sequences  generated  synthetically  for  two  representative 
geoacoustic  models  of  deep  water  sediments.  The  reconstructed  profiles 
showed  excellent  agreement  with  the  model  profiles  even  though  the  Impulse 
responses  were  lowpass  limited  to  64  Hz.  Although  the  results  of  the 
simulation  are  encouraging,  two  aspects  of  the  Inversion  method  require 
further  study:  (1)  the  effect  of  sediment  absorption  and  (2)  the  effect  of 
noise  In  the  reflection  data.  This  analysis  Is  underway,  and  a  report  on 
the  results  Is  planned  for  the  near  future. 
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